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1 Abstract 

, Linear systems with many degrees of freedom containing multiplicative and 

additive noise are considered. The steady state probability distribution for 

> 

■ equations of this kind is examined. With multiplicative white noise it is 

o ' 

shown that under some symmetry conditions, the probability distribution of 

a single component has power law tails, with the exponent independent of the 
CO , 

strength of additive noise, but dependent on the strength of the multiplicative 
J>y noise. The classification of these systems into two regimes appears to be 

^ , possible in the same manner as with just one degree of freedom. A physical 

f3 | system, that of a turbulent fluid undergoing a chemical reaction is predicted 



to show a transition from exponential to power law tails, as the reaction rate 

^ , is increased. A variety of systems are studied numerically. A replication 

H ] 

algorithm is used to obtain the Lyapunov exponents for high moments, which 
would be inaccessible by more conventional approaches. 
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I. INTRODUCTION 



This paper analyzes coupled linear equations containing additive and multiplicative noise. 
They can be written in the form 

= M ■ + A ■ + ff (1) 

= 02, • • • j 4>n) represents the variables of interest M(t) is a matrix that varies 

randomly with time. A is a time independent matrix and ff is a random time dependent 
vector. To motivate the study of eqn. fll]) we will illustrate a number of physical problems 
that are in this category. 

The motion of dye, or a temperature field, in a randomly stirred fluid leads to an equation 
for the dye density where the random velocity v is coupled to dye density in a multiplicative 
way v • V0. obeying the equation 

d t (f) + v ■ V0 = dV 2 (2) 

In two and higher dimensions this is normally studied with the incompressibility condition 
V • v = 0. d is a diffusion constant, and v = v(r,t) is a random function of position and 
time. This problem has received much recent attention [|T]-[TT|] . 

The Schroedinger equation with a random time dependent potential is another exam- 
ple [0,0, although it will not be investigated here, and this is almost identical to the 



equation for wave propagation in random media |T^|,|T3. Population growth models which 



are relevant to chemical reactions [Tj|, population biology [T_7[ and low temperature quantum 
systems are often of this form 

d t( j) = a0 + /0 + c/V 2 (3) 

where / is a random function of position and time, a and d are parameters. Dye laser theory 
also leads to an equation with both types of noise [20]. Polymers in turbulent flow also make 
use of such equations • 

The inclusion of additive noise most frequently corresponds to adding thermal fluctua- 
tions to the problem. In some cases it is a source term and replaces the effect of a boundary. 
For example in turbulent Raleigh- Bernard convection, where a liquid is heated from below, 



hot "plumes" |23j form and rise up through a liquid. In this case in eqn. (y) represents 
the temperature field, and a source term in this equation should be added representing heat 
flow through the boundaries. 

A one component problem similar to that studied here has been analyzed by Drum- 
mond [23]. However his analysis is not completely applicable to the systems mentioned 



above as in these cases the additive and multiplicative noise are uncorrelated with each 
other. However he also finds two regimes, one with all moments defined and the other with 
a divergence at a finite moment. 

The classification into different two regimes is very similar to the one component case |25[ 



where it was shown that there there are two different types of behavior found for the prob- 
ability distribution function (PDF). There it was shown that the PDF has power law tails 



in one regime and stretched exponential in another. In section pj the n-component vector 
version of these equations is considered for the case where both kinds of noise are white and 
Gaussian. Here under certain symmetry conditions it is shown that one also finds a power 
law tail for the PDF of any of the vector's components. 

An important example of non-power law tails involving many components is eqn. (^j). 
L(q) in this case should lead to a stretched exponential or an exponential tail for the PDF 
of dye density. Such distributions have been observed experimentally and have already been 
the subject of much theoretical work. The argument presented here for these tails is much 
general and just relies on mass conservation, that fact that dye is conserved. 

Is there a physical system similar to eqn. (0) where the dye particles are not conserved? 
An example of such a system would be a reactive fluid, where heat is generated by particles 
reacting with each other. Such a system would have Lyapunov exponents that cross over 
to the first regime. This leads to the interesting prediction that such systems can exhibit 
power tails. This will be discussed in section |V[ Numerical confirmation of this prediction 
is given in section |V|. 

II. DISTRIBUTION IN THE ABSENCE OF ADDITIVE NOISE 

Here we will consider eqn. ([[]) when the additive noise term rf is zero. In this case there 
is not a well defined steady state probability distribution. However this case can also be 
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understood by considering the time dependence of the distribution for long times. We shall 
see that there is a scaling form of this distribution for long times which is identical in form 
to the one-component case. 



A. Scaling of the probability distribution 

Without the additive noise, the formal solution to eqn. (|J) is 

f(t) = Te/o M W +Ad >(0) (4) 

Here T denotes the time ordered product. If this integration is discretized into time intervals 
At, then this is equivalent to the problem of the multiplication of t/ At random matrices 
of the form exp(At(Mi + A)) The matrices Mj = M(iAt) are random and independent 
because we are considering M(t) to have a white noise spectrum. 

The multiplication of random matrices has been well studied. For long times it has been 



proved that \(f>{t)\ has an exponential dependence on time for long times pqj . Fluctuations 
in this system are large in the sense that higher moments do not scale in a trivial way with 
lower moments. In one considers a single component of <fi, (fix then if 0i ~ exp(7it) , then 
<p\ ~ exp(72t), where 71 and 72 are not trivially related. This can be seen to be true even in 
the scalar case where the matrices are dimensions lxl and hence commute. In summary, 



one can define a generalized Lyapunov exponent L(q) by [p7| , p8| 

m oc e L ^ (5) 

where the brackets denote an ensemble average over the noise M. 

Now we would like to work out the form of the PDF of 0, one of the components of 0. 
This situation is similar to the problem of multifractals and its solution works out the same 
way. It can be checked that the distribution giving such scaling is 

lnP(ln0)oct/((ln0)/t) + O(^) (6) 



Another way of deriving this scaling form is by a thermodynamic analogy [p7, pq] . The 



discretized version of eqn. (HI) can be thought of as a coupled one dimensional spin system, 
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where each spin has N states but with random interactions Mi. Eqn. (|4]) is then the 
partition function of the system for a particular realization of the Mi. 

From eqn. (Q) it is clear that P(fa) does not tend towards a time independent, that 
is steady state, distribution. ((ln</>) 2 ) oc t for large t, so the the distribution continues to 
increase for small and large \fa. 



III. PROPERTIES FOR MANY COMPONENTS 



A. Equation for moments 

We will start by considering the coupled differential eqns. ([I]) where M is a random 
gaussian matrix with (Mjj(t)M^(t')) = 2T i j} t i5(t—t'), and where A is taken to be symmetric, 
real, and time independent. As will be seen below, the elements of A must be sufficiently 
negative so that there is a well defined stationary solution, that is the solution does not 
diverge for long times, ff is additive random noise that is taken to be white and uncorrelated 
(T]i(t)r]j(t')) = D5(t — t')5ij. Later in this section we will see that when M(t) is either real 
symmetric or anti-symmetric, the steady state probability distribution for fa has a power 
law tail. 



Using standard methods p9[ it is possible to derive the steady state probability distri- 
bution for cf). 

Next consider nth moments of this equation. Multiply through this equation by 
( f ) a 1 4>a 2 ■ ■ ■ 4>on where the a^s can take on any integral value between 1 and N. Integrat- 
ing this over all the fa's gives an equation relating moments of order n to moments of order 
n — 2. It is convenient to write such an equation in matrix form. Define a = (a\, «2, . . . , a n ) 
and (3 = (Pi, fa, ■ ■ ■ ,Pn) where the a^s and /3j's can take on any integral value between 1 
and N. Define <f)$ = (0 Ql a2 . . . 4> an )- Then the equation for the nth moment can be written 
in the form 

E(^-G^ = N a -?jg (8) 



where 



x 8,0 



(9) 



j 



n N n n n 

G a,0 = r a s 77ft II + <*i&<*3 Pi + T^ft^ft] II <Wa> (10) 

i=l 7=1 3 i,3 k 

j^i i<j k^i,j 

and 

n 

N(% = D $a i ,aj4 > oi\,...ai--i.,ct i+ i,...a.j- 1 ,a.j + \,...a N (H) 

i<j 

This involves moments of order n — 2. 

Eqn. (H) provides a convenient framework to analyze the general properties of moments. 
It is of the form of a matrix equation of order N n . Using this, it will be shown that these 
moments diverge at sufficiently high order if the symmetry condition, mentioned above, for 
M holds and therefore the PDF of the 0's has a power law tail. 

When there is no additive noise, that is N = 0, the Lyapunov exponent L(n) is given 
by the largest eigenvalue of the matrix G — a. 



B. Divergence of moments 

With finite additive noise in steady state, Eqn. (Q) can be rewritten 

= -£(G-a)~W N £ ( 12 ) 



The objective here is to show that there exists some finite n beyond which the moments do 
not exist. Generically this will occur when the highest eigenvalue of G — a passes through 
zero. Positive eigenvalues show that the solution has become ill defined, and that the 
moments have ceased to exist. Note that if the nth moments exist, the moments of order 
n — 2 must also exist, and hence the right hand side of eqn. QT21 ) exists. Therefore one 
accomplishes the above objective by showing that the inverse matrix on the right hand side 
of eqn. (|T2|) has eigenvalues that must pass from negative to positive, as n is increased. 
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From the end of section [III A| , the value of n where this occurs corresponds to the point 



where L(n) passes through zero. This condition is in agreement with the result obtained 
through the heuristic argument of the preceding publication |25 . 



Thus to show this divergence, one must find the dependence of the eigenvalues of G — a 
on n, the order of the moments. One must show that at least one of the eigenvalues of this 
matrix become positive for sufficiently large n. 

Note that the eigenvalues of G — a will span some finite domain E min to E max both of 
which are difficult to ascertain. However one can find a number between these two extremes 
by forming the scalar product x- (G — a) -x, where x is any unit vector. If this scalar product 
becomes positive, then at least one eigenvalue is greater than zero, and therefore the nth 
moment is divergent. If one makes the choice xoc(l,l,...,l) then the scalar product can 
be evaluated. At this point we restrict the random matrices M to be either real symmetric 
or anti-symmetric. 

22 x <x( G a,0~ a a,0) x P = 

3,0 

N N N ( _ 1 \ JV 

I E A,* + £ £<(£ M aa f) + 2^— £ M aP f) (13) 

7<5=1 7=1 a=l a/3=l 

Now we need to analyze the behavior of this expression as a function of n. The first two 
summations have prefactors proportional to n, where as the last term is proportional to 
n(n — 1) and must be positive. Therefore this scalar product must become positive for 
sufficiently large n. 

If we consider this scalar product as the Ay's are varied one obtains a similar result. 
Suppose we start with A^-'s that are sufficiently negative that the nth moment exists. As 
the Ay's become less negative, at some point the scalar product becomes positive. Therefore 
at this point the moments 05 are ill defined. Therefore for a fixed order of moments n, the 
moments can be made ill defined by varying the matrix A. 

These arguments show that for a system described by eqn (JI|), there must always be 
moments of the 0's that diverge if sufficiently high moments are examined. All of the 
moments of a given order are expected to diverge at the same time, if the matrix T is not 
transformable to block diagonal form. Block diagonal form should only occur when the 
system of eqns. (0), can be decoupled into completely separate subsystems. Aside from 



this case, all elements of (G — a) - will become divergent at the same n. This means for 
example, that (0™) should diverge at the same n as (02 _3 0l)- This argument shows that the 
divergence in moments should take place somewhere between the nth and (n+2)th moment. 
It is quite reasonable to assume that the power law tail for the PDF will take on a non- 
integral value. Note that the divergence is not dependent on the matrix N, which is a 
function of D the strength of the additive noise. Therefore the exponent of the power law 
should also be unaffected by the strength of the additive noise. However the steady state 
distribution becomes ill defined in the limit of zero additive noise as we saw in section [IT]. 
The power law tail should only depend on T and A. 

IV. PHYSICAL EXAMPLE 

Here the results of the previous section are used to understand and make predictions 
about an interesting physical system, the convection of a passive scalar field such as temper- 
ature in a random velocity flow. We will add the extra ingredient that the fluid is undergoing 
an irreversible chemical reaction A — > B. All molecules start off as type A and the reac- 
tion is exothermic. In general one expects the rate constant to be temperature dependent. 
Over a limited temperature range this dependence can be approximated as linear. It is 
easily seennicolis that this reaction adds a term acj) to the right hand side of eqn. (|2]). The 
equation under consideration is 

d t (j) + v- W(p = dW 2 <p + a<p + ri(x,t) (14) 

The last term is an additive random noise term, the correlation function in fourier space 
(|r/(k, t)\ 2 ) must go to zero as k — > to ensure heat conservation in the absence of chemical 
reactions. 

It will now be argued that this system exhibits a cross-over from an exponential tail in 
the PDF, P(T), to a power law as a is increased. 

First consider the case of no reaction, a = and no additive noise. In this case we 
will now see that L(q) must asymptote to a constant < 0. There are three ingredients in 
this argument. First the Laplacian acts as a short distance cutoff precluding variation of 
the temperature field on a smaller scale. To estimate this length scale we are assuming 
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that the velocity field has a small distance cutoff. In the case of turbulence this would be 
Kolmogorov's smallest eddy scale . Another scale of relevance is obtained by comparing the 
first two terms on the right hand side of eqn. fll4]). This yields a scale l c = d/v rms . As is 
usual with competition between convection and diffusion, it is the latter that will dominate 
below and the former that will dominate above the scale l c [30|. Therefore at a sufficiently 
small length scale the variation of T should be smooth. Second, in a closed system or one 
with zero net heat flux into the system, the integrated temperature is conserved. Third if 
the the scalar field is initially positive everywhere, it must remain so under evolution of eqn. 
(§). These three observations, of (a) conserved heat flux, (b) minimum length scale, and (c) 
positivity of 0, imply that in steady state there will be a maximum that the temperature 
can take. This would be when the entire non-zero, or non-average, temperature field is 
concentrated in a spike of size the minimum length scale. Because there is a maximum to 
<ft, it cannot grow exponentially and therefore L(q) < 0. 

But in fact from common experience one knows that without additive noise temperature 
variations will decrease in time. Physically it is apparent, that after some time an isolated 
system should come to equilibrium at a constant temperature. The form of the decrease 
in r.m.s. temperature fluctuations depends on the form of the random velocity field. The 
decrease in fluctuations will depend, for long times, on the boundary conditions. If there is 
no heat flux in or out of the system, the integrated temperature is conserved. With v = 
the process is purely diffusional, and one expects an exponential decay in fluctuations of the 
temperature field. The addition of a random velocity field should not change this conclusion. 
It should have the effect of mixing the system faster and increasing the rate of decay. If the 
boundary conditions are no longer insulating so that heat flows out of the system, the rate 
of decay should decrease further. The exact dependence is not important for the conclusions 
here, but only that in a finite size system and for a variety of boundary conditions, L(q) < 
for q > 0. 

Now consider the effect of adding a term a<fr to the right hand side of eqn. (||) to make 
the fluid reactive as described above. By letting <ft — > <pexp(at), one eliminates the term 
a<j) from the equation. Therefore the term a<p modifies the the Lyapunov exponents with 
no reaction L (q) to L(q) = L (q) + aq. Because L (q) must be convex, sufficiently large a 
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must cause L (q) + aq to cross zero at finite q. 

If we now consider the steady state behavior of P(4>) with additive noise we recall that 



because L(q) < for a = 0, P{4>) for one component models |25] must be an exponential or 
stretched exponential. For sufficiently large a, one should observe a cross over to power law 
behavior. The precise cross over behavior cannot be predicted from this analysis, only that 
such a cross over must take place. The next section confirms this cross over numerically. 



V. NUMERICAL SIMULATIONS 



We will illustrate and confirm some of the results presented here by means of computer 
simulation. Consider eqn. (|j) in one dimension. In one dimension the incompressibility 
constraint cannot be enforced and still obtain an interesting problem, so this constraint is 
dropped. The equation is discretized spatially giving 

4>i + ((j) i+1 V i+1 - (piVi) = d((f) i+1 - 20; + (15) 

This equation is solved with periodic boundary conditions, and it is clear that the total <j) 
$ = J2 (pi is conserved. Since we are interested in the fluctuations it is most convenient to 
take $ = 0. The simplest "order parameter" for this system is the total standard deviation. 

2 = |> 2 (16) 

i=i 

This is to be distinguished from the 0's being evolved by eqn. ([I]). In this section we will 
denote the order parameter simply by <fi. The different v j's are taken to be independent and 
Gaussian distributed. This equation is updated by fourth order Runge Kutta. 16 lattice sites 
are used, and the distribution P((f)) was computed by averaging. Fig. [I] plots (lnP(0))/£ 
versus (ln0)/i for two different times t. According to eqn. (^|) the resulting function should 
be independent of time, up to an overall vertical shift. The two curves have been shifted 
relative to each other and lie on top of each other within the error bars confirming the 
validity of the scaling assumed in this paper. Note that the data for the longer time PDF 
does not extend down as far as for the shorter time. There is a problem measuring the tails 
of the PDF because the number of data points becomes very small. This one dimensional 
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problem may not have an L(q) that is less than zero for positive q as the positivity of fa is 
no longer guaranteed. 

We next turn to the case where additive noise is included. We shall illustrate the tran- 
sition in behavior from an exponential tail in the PDF to a power law tail in the case of a 
reactive randomly stirred fluid. Adding a term a<pi to the right hand of eqn. QT5| ) repre- 



sents a one dimensional discretized version of the reacting system described in section |TV . 
Periodic boundary conditions are still employed here but we start with the condition that 
the sum over all sites of fa is zero. In this case for a finite size system and a = 0, (fa) will 
decrease exponentially. This is seen from the fact that when the velocity field is zero, the 
system is diffusive and fa decreases exponentially. When velocity is added, this does not 
alter this result. From the previous section it was predicted that a system with sufficiently 
large a should have power law tails. 

Fig. (0) shows <f) versus P(fa) for different values of a. If a is too large all the L(q) 
are positive and the system is unstable. For a too small there is no discernible power law 
behavior. However this cannot be ruled out based on the data. There is a narrow region 
where clear power law tails are manifest. 

The method just used for obtaining the PDF has the draw back of not being able to 
probe the tails of the distribution. Equivalently, the high order moments and their L(g)'s 
are difficult to obtain. It is the convergence of high order moments that presents a problem. 
Multiplicative noise equations like eqn. (|IJ) have highly intermittent behavior. This is a result 
of the fact that different moments have different Lyapunov exponents. Therefore the typical 
behavior of a high moment will be very different from its true average. Flucutations that 
are exponentially unlikely in time dominate the correct answer. Therefore straightforward 
averaging is practically limited to very low q, before the number of number of runs that are 
needed to obtain convergence becomes far too large. 

A way around this problem is to use a "replication algorithm". Similar methods have 



been used in quantum simulations for the past fifteen years ||31|| . In order for this method to 
work, one needs the order parameter <p to be positive definite, which it is in most situations, 
as in the above example. A large number of copies of the system described by eqn. ([I]) 
are made. These are all updated to the next time step, with independently chosen random 
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numbers, the matrices M. A weight w is computed by taking the ratio w = (f) n (t + At)/(ft n (t). 
Again can be taken to be the total <fi field, or the standard deviation of it, as above. Using 
this total rather than 0; at a single point i leads to better convergence and is therefore 
what is done in practice. Once the weights w are determined, the integer part of w is taken 
and that number of copies is made. The remainder of w is taken into consideration, by 
generating a random number between and 1 and making an additional copy if the random 
number is smaller than the remainder. Additionally, to stop the exponential explosion in 
the number of copies, the weights are all normalized so that the number of copies is kept 
almost constant. 

The above method takes into account the idea of importance sampling. When one is 
interested in the average value of the moment <p n+1 (t), this is not computed directly, but 
the average of <p is computed. It is computed for a biased ensemble of configurations, 
each one having a probability proportional to <p n of appearing. Therefore one is comput- 
ing (0 n 0)/(0 n ). This method is most powerful with no additive noise. In this case the 
exponential divergence of this ratio is computed for long times, which gives L(n+ 1) — L(n). 

The method is illustrated first for the two dimensional version of eqn. (0) with no additive 
noise. To make sure that V • v = 0, v is derived from a vector potential v = V x A. The 
velocity field v is taken to be smooth but randomly varying in time. A is taken to be of the 
form 

A( X ,y,t) = M( t )sin( 2T( "7° (t)) )sin( 2 " to 7° (t)) )i (17) 

New values of x® and yo are randomly chosen every At steps. M(t) is a Gaussian random 
variable. This describes a circular motion with periodic boundary conditions, moving with 
a random angular velocity. There is no problem obtaining satisfactory convergence even for 
large moments. The value of In as a function of time is plotted for three different values 
of q in fig. ^ for an 8 x 8 system.. Approximately 720 copies are run in this simulation. 
The slopes of these lines give the difference in Lyapunov exponents AL = L(q + 1) — L(q). 
This difference is plotted versus q in fig. | in a system of size 8x8. The same quantities 
are plotted in fig. || in a system of size 16 x 16. All other parameters are kept the same 
and 180 copies were updated. Note that the variation of L(q) is much smaller for this larger 
system. Also note that L(q) is always negative for positive q and appears convex. By adding 
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a chemical reaction of the type described in section [TV] for sufficiently large a, L(q) should 
cross zero leading to power law tails. 

The final example considered here is the population growth model, eqn. (H), with a = 0. 
This has been studied extensively because of its relevance to quantum mechanical many 
body systems. The Lyapunov exponent L(q) corresponds to the ground state energy of q 
particles interacting with each other via a two body potential as was first shown by Sugiyama 



and Koonin [Tj|. Here we are considering a lattice model so that r takes discrete values and 
V 2 is the discrete Laplacian. Without the last term, the above equation is just the diffusion 
equation, or the imaginary time Schroedinger equation for a free particle. The inter-particle 
interaction is a result of the multiplicative term as we shall see shortly. Assume that / is 
Gaussian noise with zero mean and the correlation function (/(r, t)f(r', t')) = v(r — r')5(t — 
t>). 

For discrete r consider the nth moment averaged at equal times 

^(n,r 2 ,...,r n ,t) = (0(n,t)0(r 2 ,t)...0(r n ,t)). (18) 



In the notation of section p| , ^ = (fi$ and a corresponds to the coordinates (r 1; r 2 , . . . , r n ) 
Translating eqns. (|8|-p!0D into these new symbols it can be seen that \1/ obeys the many body 
Schroedinger equation with imaginary time. 

-4? = (19) 

H = K + U is the Hamiltonian of the system which is the sum of the kinetic and potential 
energy, K and U respectively 

n 

K = -J]V 2 (20) 



t=i 

and 

71 n 
U = - 5>( ri - r,) - -«(0) (21) 

i<3 Z 

Note that u need not be positive as / may be chosen to be complex. This allows for the 
simulation of systems of particles with both attractive and repulsive interactions. Also 
note that the last term in the potential energy is constant and therefore does not present 
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any difficulties. The addition of a term — U ex t(r)<j)(r) to eqn. (|3]) incorporates an external 
potential U ext in the Hamiltonian H. 

Thus the wavefunction of a system of n bosons, with certain initial conditions, is simply 
related to the nth moment of equation (|3]). To obtain the ground state energy of n bosons, 
E n , one need calculate the Lyapunov exponent L(n) = —E n . One is normally not interested 
in the constant term — f'u(O) that appears in the expression for the potential energy, so this 
is subtracted out. Eqn. (|3|) has been studied numerically by path integral techniques [JTSj . 
It is worth noting that the replication algorithm presented here is also an efficient procedure 
for obtaining ground state energies. Again, straightforward simulation of this equation will 
not work because the higher order moments will not converge in a reasonable time. The 
replication algorithm efficiently solves this problem. 

The method with the addition of the replication algorithm can handle large number 
of particles and is quite computationally efficient. Fig. ^| shows the difference in energies 
fi n = E n — E n+ i as a function of n for a 10 x 10 system with an on-site attractive potential 
U = —0.5. Approximately 900 copies are run in this simulation. In two dimensions there 
is a transition between localized and delocalized particles as a function of density. The two 
dimensional transition can be understood as a competition between kinetic and potential 
energy. The kinetic energy scales as n/R 2 where R is the correlation length. The potential 
energy scales as un(n — l) 2 / (2ttR 2 ). In the delocalized phase R is the system size L, and in 
the localized phase R should become as small as possible which is one lattice spacing. This 
corresponds to all particles occupying the same lattice site. The transition between localized 
and delocalized phases will then occur at approximately n = 2tc/u. For the figure shown 
the transition occurs within a factor of two from the simple estimate given above. This 
transition is not a thermodynamic transition as it occurs at a constant number of particles, 
almost independent of volume. 

An interesting feature of the replication algorithm applied to this system is that it exhibits 
metastability and hysteresis right above the transition. If the initial field 0(r) is delocalized 
then it will stay so for many iterations of the replication algorithm. It will eventually jump 
to the localized phase and remain in that state. The interpretation of this is similar to 
nucleation theory ]32[ . 
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VI. CONCLUSIONS 



The PDF for eqn. (|1]) has been analyzed for many components. It was found that PDF 
will exhibit power law tails when a symmetry condition holds for the random multiplicative 
matrix M. 

A passive scalar such as temperature in a random velocity field is predicted to have a PDF 
with an exponential or stretched exponential tail in agreement with earlier experimental and 
theoretical work PJ3H1TI . If the fluid is reactive this is predicted to become power law when 
the reaction rate is sufficiently large. Numerical evidence was presented supporting this 
prediction. The numerical work suggests that this power law behavior is most easily seen in 
small systems where the velocity correlation length is comparable to system size. 

It is interesting that the multi-component Gaussian case should exhibit the same behavior 
as the a one component system albeit one that is no longer gaussian. An argument can be 
given to understand this equivalence with no additive noise. A multidimensional system has 
a variety of relaxation times. Consider the longest of these T re \. For times much longer than 
T re i a system looks statistically identical to how it did originally except for an overall shift in 
scale. To make this concrete, suppose we shift the scale periodically after a time AT >> T re i 
so that \<j)\ 2 = 1. The shift in scale will not be completely deterministic but will fluctuate 
because of the stochastic nature of the equation. For a long time MAT, the overall scale 
factor is obtained by multiplying the M shifts in scale together. Each of these scale shifts 
is random and in general is expected to be non-gaussian. The nature of this non-gaussian 
distribution depends on the evolution of the system within a relaxation time. Therefore the 
long time evolution of this system can still be properly modeled by a one component system 
as was investigated previously p5 
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FIGURES 

FIG. 1. The probability distribution of a field <fi obeying eqn. (15) in the absence of additive 
noise. The system is in one dimension and is plotted with rescaled axis. The horizontal axis is 
\a.<j)/t and the vertical axis is (lnP)/i. Two different times t\ = 12.8 and t<i = 38.4 are shown. The 
variance of the multiplicative noise is 2.0, the diffusion coefficient d is 0.1. 

FIG. 2. The probability distribution of field <fi obeying eqn. (14) in one dimension, Six different 
values of a are shown. The most steeply descending curve is for a = 0, and in order order of 
decreasing steepness, a = 0.40, a = 0.47, a = 0.49, a = 0.50, with the top curve at a = 0.51. 
The equation was discretized and the number of lattice sites was chosen to be 8. The variance of 
the both the multiplicative and additive noise is 0.6. 

FIG. 3. The logarithm of {(ft) as a function of time, for a two dimensional passive scalr field. 
The average is a weighted average corresponding to the weights cf) 1 and 4> s and calculated 
using the replication algorithm described in the text. The slopes of these lines correspond to 
AL = L(q + 1) — L(q), with q = 1, 8 and 16. The bottom curve corresponds to q = 1, the middle 
curve to q = 8, and the top curve is for q = 16 

FIG. 4. The difference in Lyapunov exponent AL as a function of q for the two dimensional 
motion of a passive scalar field on an 8 x 8 lattice. Note all of the values are monotonically 
decreasing confirming that L(q) is convex. 

FIG. 5. The difference in Lyapunov exponent AL as a function of q for the same system as the 
previous figure but on a 16 x 16 lattice. Note the variation of L(q) with q is much smaller in this 
larger system. 

FIG. 6. The chemical potential at zero temperature, ^ = E n — E n+ \ as a function of n, the 
number of bosons in a two dimensional 10 x 10 lattice. The on-site attractive energy is —0.5. The 
energies are calculated by the replication method described in the text. Note the discontinuity at 
n = 19. 
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Fig. 5 
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Fig. 6 
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